# Clears work environment
rm(list = ls())

# Sets working directory
setwd('C:/Users/Jason/Box Sync/Home Folder jdt34/733 - Maximum Likelihood Estimation/Lai & Slater 2006 replication')

# Loads data and libraries
load(file='cleanedLSdata.RData')
loadPkg(packs)

# Creates coefficient plot for Lai & Slater model with fundamental uncertainty
betas.nb = data.frame(t(apply(negbdraws, 2, FUN=cimaker)))
rownames(betas.nb) = NULL
colnames(betas.nb) = c('low95', 'low90', 'betaest', 'upp90', 'upp95')
betas.nb$parameter = factor(c('(Intercept)', IVs), levels=rev(c('(Intercept)', IVs)))
betas.nb$sig = 0
betas.nb$sig[betas.nb$low95>0 | betas.nb$upp95<0] = 1
ggcoefnb = ggplot(data=betas.nb, aes(x=parameter, y=betaest, ymin=low95, ymax=upp95, color=factor(sig))) + 
  geom_hline(y=0, linetype=2, color='darkgrey') +
  geom_errorbar(width=.25) + 
  geom_linerange(size=1, aes(ymin=low90, ymax=upp90)) +
  geom_point() +
  guides(color=F) +
  xlab('') +
  ylab('') +
  coord_flip()

# Creates coefficient plot for zero-inflated model with fundamental uncertainty
betas.zi = data.frame(t(apply(zinbdraws, 2, FUN=cimaker)))
rownames(betas.zi) = NULL
colnames(betas.zi) = c('low95', 'low90', 'betaest', 'upp90', 'upp95')
betas.zi$parameter = factor(c('(Intercept)', IVs), levels=rev(c('(Intercept)', IVs)))
betas.zi$model = factor(c(rep('Count process', 10), rep('Zero process', 10)), levels=c('Count process', 'Zero process'))
betas.zi$sig = 0
betas.zi$sig[betas.zi$low95>0 | betas.zi$upp95<0] = 1
ggcoefzi = ggplot(data=betas.zi, aes(x=parameter, y=betaest, ymin=low95, ymax=upp95, color=factor(sig))) + 
  geom_hline(y=0, linetype=2, color='darkgrey') +
  geom_errorbar(width=.25) + 
  geom_linerange(size=1, aes(ymin=low90, ymax=upp90)) +
  geom_point() +
  scale_y_continuous(breaks=c(-15, -10, -5, 0, 5)) +
  facet_wrap(~model, ncol=2) +
  guides(color=F) +
  xlab('') +
  ylab('') +
  coord_flip(ylim=c(-17.5, 7.5))

ggsave(file='nbcoefs.pdf', ggcoefnb, width=9, height=3, units='in')
ggsave(file='zicoefs.pdf', ggcoefzi, width=9, height=3, units='in')

save(file='panel2_coefficient_plots.RData', betas.nb, ggcoefnb, betas.zi, ggcoefzi)
